Raster parallel process in R

Install libraries if there are no already in your computer.

## [1] "raster"

## Loading required package: raster

## Warning: package 'raster' was built under R version 4.0.5

## Loading required package: sp

## [1] "terra"

## Loading required package: terra

## Warning: package 'terra' was built under R version 4.0.5

## terra version 1.1.17

## [1] "snow"

## Loading required package: snow

## Warning: package 'snow' was built under R version 4.0.4

## [1] "parallel"

## Loading required package: parallel

## 
## Attaching package: 'parallel'

## The following objects are masked from 'package:snow':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, clusterSplit, makeCluster, parApply,
##     parCapply, parLapply, parRapply, parSapply, splitIndices,
##     stopCluster

## [1] "foreach"

## Loading required package: foreach

## Warning: package 'foreach' was built under R version 4.0.5

## [1] "rasterVis"

## Loading required package: rasterVis

## Warning: package 'rasterVis' was built under R version 4.0.5

## Loading required package: lattice

## Loading required package: latticeExtra

## Warning: package 'latticeExtra' was built under R version 4.0.5

## [1] "mapview"

## Loading required package: mapview

## Warning: package 'mapview' was built under R version 4.0.4

## GDAL version >= 3.1.0 | setting mapviewOptions(fgb = TRUE)

## [1] "ggplot2"

## Loading required package: ggplot2

## 
## Attaching package: 'ggplot2'

## The following object is masked from 'package:latticeExtra':
## 
##     layer

## [1] "lubridate"

## Loading required package: lubridate

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:terra':
## 
##     intersect, union

## The following objects are masked from 'package:raster':
## 
##     intersect, union

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

## [1] "pryr"

## Loading required package: pryr

## Warning: package 'pryr' was built under R version 4.0.5

## Registered S3 method overwritten by 'pryr':
##   method      from
##   print.bytes Rcpp

## 
## Attaching package: 'pryr'

## The following object is masked from 'package:terra':
## 
##     dots

## The following object is masked from 'package:raster':
## 
##     subs

Download files, alternative is possible use the function raster::getData. More details on https://worldclim.org/data/worldclim21.html

#download worldclim climate data
download.file('https://biogeo.ucdavis.edu/data/worldclim/v2.1/base/wc2.1_10m_tavg.zip',temp, mode="wb")
#raster::getData()
unzip(temp, exdir=tempd)
#dir(tempd) # explore the temp directory

Explore the data

file_names = list.files(tempd,'wc2.1_10m_tavg')
tmean <- stack(file.path(tempd, file_names))
mapview(tmean[[7]])

## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 2332800 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  2332800 '

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition

## Warning in rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, : 104
## projected point(s) not finite

## Warning in rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, : 104
## projected point(s) not finite

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition

## Paraleling Create cluster and load functions. Alternative you can create your own scrip file with function and just load this file.

cl <- makeCluster(4, type = "SOCK")
out=clusterEvalQ(cl,library(raster))
clusterExport(cl, list('my_function'))
#out =clusterEvalQ(cl, source('my_functions.R'))

We going to use two way of paralleling in R with the snow/parallel and foreach libraries.

system.time({
  clusterExport(cl,list('tmean'))
  mean_planet = parLapply(cl,c(1:12), function(l){
    my_function(tmean[[l]])
  })
  mean_planet = unlist(mean_planet)
})

##    user  system elapsed 
##    0.02    0.01    2.05

system.time({
  mean_planet2 = foreach(l = c(1:12),
                         .combine = rbind,
                         .errorhandling = 'remove') %dopar%
    my_function(tmean[[l]])
}) 

## Warning: executing %dopar% sequentially: no parallel backend registered

##    user  system elapsed 
##    3.48    0.53    4.02

results = data.frame(layer_name = names(tmean),
                     month =  month(1:12,label=TRUE, abbr = TRUE),
                     parLapply_out =mean_planet,
                     foreach_out = mean_planet2)

The first plot show the global mean temperature for each month. The second figure show that the two paralleling outcome the same results

ggplot(results,aes(x=month,y=parLapply_out, col='parLapply_out'))+
  geom_point()+labs(y='Mean temperature (°C)', col='', x='')+
  geom_point(aes(y=foreach_out,col = 'foreach_out'))+theme_minimal()

Split large raster

We download the raster and make a split overlapped extent. The overlapped extent is only for the use of focus function, any other process should use a non-overlapped extent! Also is possible use the function focus in raster. but it is slower function.

file_raster = 'https://www.dropbox.com/s/onp3f75pab7z3j9/intro-1508529851.jpg?dl=1'

r_raster = stack(file_raster)
e =extent(r_raster)
extent_r =list(r1 =extent(0, floor(e@xmax/2)+2, 0, floor(e@ymax/2)+2),
               r2 = extent(ceiling (e@xmax/2), e@xmax, floor (e@ymax/2), e@ymax),
               r3 =  extent(0, floor (e@xmax/2)+2,  floor (e@ymax/2), e@ymax),
               r4 =  extent(ceiling (e@xmax/2), e@xmax,  0, floor (e@ymax/2)+2))

We call the parLapply from a lapply environment, so we need export elements from this environment (env).

tmean_smooth = lapply(c(1:3) , function(l){
  env= where('l')
  clusterExport(cl, list('extent_r','r_raster','l'),envir=env)
  r_parts = parLapply(cl,seq_along(extent_r),function(part){
    r = r_raster[[l]]
    r_part <- crop(r, extent_r[[part]])
    #focal(r_part,w=matrix(1,3,3), fun='mean',na.rm=TRUE)
    terra::focal(r_part,w=matrix(1,3,3), fun='mean',na.rm=TRUE)
    #r_part
  })
  
  m1 <- mosaic(r_parts[[1]], r_parts[[2]],r_parts[[3]],r_parts[[4]], fun=mean)
  
  m1
})
#r_smooth = stack(tmean_smooth[[1]],tmean_smooth[[2]],tmean_smooth[[3]])
r_smooth = do.call(stack,tmean_smooth)

plotRGB(r_smooth)

stopCluster(cl)

Versus original

plotRGB(r_raster)